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Abstract 

A  new  high  order  finite-difference  method  utilizing  the  idea  of  Harten  ENO 
subcell  resolution  method  is  proposed  for  chemical  reactive  flows  and  combustion. 

In  reaction  problems,  when  the  reaction  time  scale  is  very  small,  e.g.,  orders  of 
magnitude  smaller  than  the  fluid  dynamics  time  scales,  the  governing  equations 
will  become  very  stiff.  Wrong  propagation  speed  of  discontinuity  may  occur  due 
to  the  underresolved  numerical  solution  in  both  space  and  time.  The  present 
proposed  method  is  a  modified  fractional  step  method  which  solves  the  convection 
step  and  reaction  step  separately.  In  the  convection  step,  any  high  order  shock¬ 
capturing  method  can  be  used.  In  the  reaction  step,  an  ODE  solver  is  applied 
but  with  the  computed  flow  variables  in  the  shock  region  modified  by  the  Harten 
subcell  resolution  idea.  For  numerical  experiments,  a  fifth-order  finite-difference 
WENO  scheme  and  its  anti-diffusion  WENO  variant  are  considered.  A  wide  range 
of  ID  and  2D  scalar  and  Euler  system  test  cases  are  considered.  Studies  indicate 
that  for  the  considered  test  cases,  the  new  method  maintains  high  order  accuracy 
in  space  for  smooth  flows,  and  for  stiff  source  terms  with  discontinuities,  it  can 
capture  the  correct  propagation  speed  of  discontinuities  in  very  coarse  meshes  with 
reasonable  CFL  numbers. 

Key  words:  stiff  reaction  term,  shock  capturing,  detonation,  WENO,  ENO  subcell 
resolution 

*  Department  of  Mathematics  and  Statistics,  Florida  International  University,  Miami,  FL  33199 
iDivision  of  Applied  Mathematics,  Brown  University,  Providence,  RI  02912 
1NASA  Ames  Research  Center,  Moffett  Field,  CA  94035 
^Lawrence  Livermore  National  Laboratory,  Livermore,  CA  94551 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
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1  Introduction 


In  simulating  hyperbolic  conservation  laws  in  conjunction  with  an  inhomogeneous  stiff 
source  term,  if  the  solution  is  discontinuous,  spurious  numerical  results  may  be  produced 
due  to  different  time  scales  of  the  transport  part  and  the  source  term.  This  numerical 
issue  often  arises  in  combustion  and  high  speed  chemical  reacting  flows. 

The  reactive  Euler  equations  in  two  dimensions  have  the  form 

Ut  +  F(U)x  +  G(U)y  =  S(U),  (1) 

where  U,  F(U ),  G(U)  and  S(U)  are  vectors.  If  the  time  scale  of  the  ordinary  differential 
equation  (ODE)  Ut  =  S(U)  for  the  source  term  is  orders  of  magnitude  smaller  than 
the  time  scale  of  the  homogeneous  conservation  law  Ut  +  F(U)X  +  G{U)y  =  0  then 
the  problem  is  said  to  be  stiff.  In  high  speed  chemical  reacting  flows,  the  source  term 
represents  the  chemical  reactions  which  may  be  much  faster  than  the  gas  flow.  This  leads 
to  problems  of  numerical  stiffness.  Insufficient  spatial/temporal  resolution  may  cause 
an  incorrect  propagation  speed  of  discontinuities  and  nonphysical  states  for  standard 
dissipative  numerical  methods. 

This  numerical  phenomenon  was  first  observed  by  Colclla  et  al.  [13]  in  1986  who  con¬ 
sidered  both  the  reactive  Euler  equations  and  a  simplified  system  obtained  by  coupling 
the  inviscid  Burgers  equation  with  a  single  convection/reaction  equation.  LeVeque  and 
Yee  [23]  showed  that  a  similar  spurious  propagation  phenomenon  can  be  observed  even 
with  scalar  equations,  by  properly  defining  a  model  problem  with  a  stiff  source  term. 
They  introduced  and  studied  the  simple  one-dimensional  scalar  conservation  law  with 
an  added  nonhomogeneous  parameter  dependent  source  term 


ut  +  ux  =  S(u), 

(2) 

S(u )  =  -fm{u  -  i)(w  -  1), 

(3) 

where  the  parameter  1  can  be  described  as  the  reaction  time.  When  /i  is  very  large,  a 
wrong  shock  speed  phenomenon  will  be  observed  in  a  coarse  mesh.  In  order  to  isolate 
the  problem,  LeVeque  and  Yee  solve  (2)  and  (3)  by  the  fractional  step  method.  For  the 
particular  source  term,  the  reaction  (ODE)  step  of  the  fractional  step  method  can  be 
solved  exactly.  They  found  that  the  propagation  error  is  due  to  the  numerical  dissipation 
contained  in  the  scheme,  which  smears  the  discontinuity  front  and  activates  the  source 
term  in  a  nonphysical  manner.  By  increasing  the  spatial  resolution  by  an  order  of 
magnitude,  they  were  able  to  improve  the  correct  propagation  speed. 

It  is  noted  that,  in  a  general  stiff  source  term  problem,  a  sufficient  spatial  resolution 
is  as  important  as  temporal  resolution  when  the  reaction  step  of  the  fractional  step 
method  cannot  be  solved  exactly.  A  study  linking  spurious  numerical  standing  waves 
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for  (2)  and  (3)  by  first  and  second-order  spatial  and  temporal  discretizations  can  be 
found  in  Lafon  and  Yee  [22,  21]  and  Griffiths,  Stuart  and  Yee  [15]. 

For  the  last  two  decades,  this  spurious  numerics  has  attracted  a  large  volume  of 
research  work  in  the  literature  (see,  e.g.,  [5,  27,  6,  31,  8,  15,  24,  1,  7,  26])  Various 
strategies  have  been  proposed  to  overcome  this  difficulty.  Since  numerical  dissipation 
that  spreads  the  discontinuity  front  is  the  cause  of  the  wrong  propagation  speed  of 
discontinuities,  a  natural  strategy  is  to  avoid  any  numerical  dissipation  in  the  scheme. 
In  combustion,  level  set  and  front  tracking  methods  were  used  to  track  the  wave  front  to 
minimize  this  spurious  behavior  [24,  1,  7,  26].  In  [11,  12],  Chorin  introduced  the  random 
choice  method  which  is  based  on  the  exact  solution  of  Riemann  problems  at  randomly 
chosen  locations  within  the  computational  cells  and  does  not  need  to  introduce  any 
viscosity.  It  has  been  successfully  used  in  [13,  25]  for  the  solution  of  underresolved 
detonation  waves.  However,  it  is  difficult  to  eliminate  all  numerical  viscosity  in  a  shock- 
capturing  scheme.  There  are  also  many  works  on  modifying  shock-capturing  methods 
for  this  problem  in  the  literature.  Fractional  step  methods  are  commonly  used  for 
allowing  an  underresolved  meshsize.  Such  methods  solve  the  homogeneous  conservation 
law  (i.e.,  the  convection  step)  and  the  ODE  system  (i.e.,  the  reaction  step)  separately. 
In  [9,  10],  Chang  applied  Harten’s  subccll  resolution  method  [16]  in  a  finite  volume  ENO 
method  in  the  convection  step  with  exact  time  evolution,  which  is  able  to  produce  a  zero 
viscosity  shock  profile  in  the  nonreacting  flow.  The  time  evolution  is  advanced  along 
the  characteristic  line.  Correct  results  were  shown  in  the  one-dimensional  scalar  case. 
However  it  seems  difficult  to  extend  this  method  to  one- dimensional  systems  or  multi¬ 
dimensional  scalar  equations  or  systems,  due  to  the  requirement  of  exact  time  evolution. 
In  [14],  Engquist  and  Sjogreen  proposed  a  simple  temperature  extrapolation  method 
based  on  a  finite  difference  ENO  scheme  with  implicit  Runge-Kutta  time  discretization, 
which  uses  a  first /second  order  extrapolation  of  the  temperature  value  from  outside  the 
shock  profile.  Their  approach  is  easily  extended  to  multi-dimensions.  However,  their 
method  is  not  a  fractional  step  method,  and  it  does  not  work  well  in  the  situation  of 
insufficient  spatial  resolution.  Hezcl  et  al.  [17]  presented  a  modified  fractional  step 
method  for  detonation  waves  in  which  the  exact  Riemann  solution  is  used  to  determine 
where  burning  should  occur.  Bao  and  Jin  [2,  3,  4]  proposed  a  random  projection  method 
based  on  the  fractional  step  method  where  in  the  convection  step  a  standard  shock- 
capturing  scheme  is  used,  and  in  the  reaction  step  a  projection  is  performed  to  make 
the  ignition  temperature  random.  They  have  successfully  applied  this  method  to  various 
problems  in  one-  and  two- dimensions.  However  they  assume  an  a  priori  stiff  source.  In 
[32],  Tosatto  and  Vigevano  proposed  a  MinMax  scheme,  which  is  based  on  a  two- value 
variable  reconstruction  within  each  cell,  where  the  appropriate  maximum  and  minimum 
values  of  the  unknown  are  considered.  The  scheme  may  be  applied  with  no  difficulties  to 
both  stiff  and  nonstiff  problems.  Only  one-dimensional  problems  were  tested.  However, 
it  seems  difficult  to  generalize  either  the  random  projection  method  or  the  MinMax 
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method  to  higher  order  accuracy. 

Our  objective  in  this  study  is  to  develop  a  high  order  finite  difference  method  which 
can  capture  the  correct  detonation  speed  in  an  underresolved  mesh  and  will  maintain 
high  order  accuracy  in  the  smooth  part  of  the  flow.  The  first  step  of  the  proposed 
fractional  step  method  is  the  convection  step  which  solves  the  homogeneous  hyperbolic 
conservation  law  in  which  any  high-resolution  shock-capturing  method  can  be  used. 
The  aim  in  this  step  is  to  produce  a  sharp  wave  front,  but  some  numerical  dissipation 
is  allowed.  The  second  step  is  the  reaction  step  where  an  ODE  solver  is  applied  with 
modified  transition  points.  Here,  by  transition  points,  we  refer  to  the  smeared  numerical 
solution  in  the  shock  region,  which  is  due  to  the  dissipativity  of  a  shock-capturing 
scheme.  Because  the  transition  points  in  the  convection  step  will  result  in  large  erroneous 
values  of  the  source  term  if  the  source  term  is  stiff,  we  first  identify  these  points  and 
then  extrapolate  them  by  a  reconstructed  polynomial  using  the  idea  of  Harten’s  subcell 
resolution  method.  Unlike  Chang’s  approach,  we  apply  Harten’s  subcell  resolution  in 
the  reaction  step.  Thus  our  approach  is  flexible  in  allowing  any  shock-capturing  scheme 
as  the  convection  operator.  In  the  reaction  step,  since  the  extrapolation  is  based  on 
the  high  order  reconstruction,  high  order  accuracy  can  be  achieved  in  space.  The  only 
drawback  in  our  current  approach  is  that  the  temporal  accuracy  will  only  be,  at  most, 
second-order  due  to  the  time  splitting,  which  is  common  for  most  of  the  previous  methods 
for  stiff  sources. 

This  paper  is  organized  as  follows.  Section  2  introduces  the  proposed  fractional  step 
method  with  subcell  resolution  for  the  one-dimensional  scalar  model  problem  in  [23]. 
The  high  order  accuracy  of  the  method  and  its  capability  of  capturing  the  correct  speed 
of  propagation  of  discontinuity  are  illustrated  with  numerical  examples.  The  proposed 
method  is  extended  to  two-dimensional  scalar  problems  with  numerical  examples  in  Sec¬ 
tion  3.  In  Section  4  the  method  is  extended  to  one- dimensional  reactive  Euler  equations 
with  a  one-step  reaction.  The  numerical  examples  include  the  Chapman-Jouguet  (C-J) 
detonation  waves.  The  method  is  extended  to  two-dimensional  reactive  Euler  equa¬ 
tions  in  Section  5  and  demonstrated  by  one  numerical  example.  Section  6  contains  the 
conclusion  and  remarks  on  future  work. 


2  Numerical  method  for  ID  scalar  problems 


We  first  introduce  the  proposed  method  for  the  scalar  model  problem  in  [23],  i.e. , 

ut  +  f(u)x  =  S(u )r 
S(u)  =  —/it/.  ( u  —  a)  (w  —  1), 


(4) 

(5) 
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with  the  initial  condition 


(6) 


u(x,  0) 


1,  x  <  x0 

0,  x  >  x0  ’ 


where  a  is  a  parameter,  0  <  a  <  1,  and  Xq  is  the  position  of  the  initial  discontinuity. 

The  general  fractional  step  approach  is  as  follows.  The  numerical  solution  at  time 
level  tn+ 1  is  approximated  by 


un+1  =  R(At)A(At)un.  (7) 

The  convection  operator  A  is  defined  to  approximate  the  solution  of  the  homogeneous 
part  of  the  problem  on  the  time  interval,  i.e., 

th  +  f(u)x  =  0,  tn  <  t  <  tn+1.  (8) 

The  reaction  operator  R  is  defined  to  approximate  the  solution  on  a  time  step  of  the 
reaction  problem: 

du  . 

—  =  S(u),  tn<t<tn+1.  (9) 

In  the  Strang-splitting  in  [30],  the  numerical  solution  at  time  step  tn+ 1  is  computed  by 

«"+1  =  A  (A)  R(At)A  (10) 

where  the  convection  operator  is  over  a  time  step  At  and  the  reaction  operator  is  over 
At/2.  The  two  half-step  reaction  operations  over  adjacent  time  steps  can  be  combined 
to  save  cost. 

Next,  we  introduce  the  proposed  fractional  step  methods  for  the  convection  step  and 
the  reaction  step  separately. 

2.1  Convection  operator 

Any  high  resolution  shock  capturing  operator  can  be  used  in  the  convection  step.  The 
purpose  in  this  step  is  to  minimize  the  transition  points  in  the  shock  region.  In  this 
paper,  we  use  the  framework  of  high  order  finite  difference  WENO  schemes  [19]  with  a 
TVD  Runge-Kutta  time  discretization  to  solve  the  one- dimensional  scalar  conservation 
law 

ut  +  f(u)x  =  0.  (11) 

In  particular,  for  the  scalar  case,  we  are  interested  in  applying  the  anti- diffusive  flux 
corrections  [33]  for  the  WENO  scheme  to  obtain  sharp  resolution  for  contact  disconti¬ 
nuities.  In  this  section,  we  will  briefly  introduce  this  anti- diffusive  WENO  scheme  for 
Eq.  (11). 
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Let  Xi,  i  —  1, . . . ,  N  be  a  uniform  (for  simplicity)  mesh  of  the  computational  domain, 
with  mesh  size  Ax  =  Xi+\  —  Xi.  An  explicit  conservative  fully  discrete  finite  difference 
scheme  has  the  form 

<+1  -  <  +  Hff+U 2  -  A 1/2)  =  0,  (12) 

where  u ”  is  the  approximation  to  the  point  value  u(xi,tn),  A  =  At/ Ax,  and  f™+1/ 2  is  the 
numerical  flux. 

2.1.1  Modified  TVD  Runge-Kutta  time  discretization  for  anti-diffusive  WENO 
schemes 

The  third-order  TVD  Runge-Kutta  (RK3)  time  discretization  [28]  can  be  written  as 


«(1)  =  un  +  AtL(un),  (13) 

u(2)  =  un  +  ^AtL(un)  +  ^AtL(uw),  (14) 

un+ 1  =  un  +  -AtL(un)  +  -AtL(u{1))  +  -AtL(ui2)),  (15) 

6  6  3 

where  L  is  the  spatial  discretization  of  —f{u)x.  The  modified  Runge-Kutta  time  dis¬ 
cretization  for  anti-diffusive  WENO  schemes  [33]  is  given  by 

w(1)  =  un  +  AtL^(un),  (16) 

w(2)  =  un+^AtL^(un)  +  ^AtL^(uw),  (17) 

un+ 1  =  un+-AtL^(un)  +  -AtLw(uw)  +  -AtLw(u{2)),  (18) 

6  6  3 

where  the  operators  L ^  are  defined  as 

=  -(Al/2-Al/2)/^.  *  =  1.2,3,  (19) 


with  the  anti-diffusive  flux  and  the  modified  anti- diffusive  fluxes  ,2  and 

defined  in  the  next  subsection. 


2.1.2  Anti-diffusive  flux  with  high  order  WENO  finite  difference  reconstruc¬ 
tion 

The  anti- diffusive  flux  for  WENO  scheme  with  RK3  is  defined  by 


7(i) 
Ji+ 1/2 


fi+ 1/2  +  0jminmod 


Uj 


H'i— 1 


+  fi- 


1/2 


f i+1/2’  f *+1/2  fi 


*+1/2 


(20) 
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where 


minmod(a,  b )  =  sgn(a)  ■  max{0,  min[|a|,  b  sgn(a)]},  (21) 

f~+ 1/2  and  f++1/ 2  are  the  two  upwind  biased  fluxes  based  on  WENO  stencils  with  one 
more  point  to  the  left  and  to  the  right,  respectively.  For  WENO-Roe  schemes,  the 
numerical  flux  is  chosen  as  f//+l/0  for  f'(u)  >  0  and  /+. 1/2  for  f(u)  <  0. 

The  function  0  is  a  discontinuity  indicator  which  is  close  to  0  in  smooth  regions  and 
close  to  1  near  a  discontinuity.  The  indicator  0  in  [33]  is 


0i 


A 

Pi  +  7  i 


(22) 


where 


0—1  O  T  O  Pi 


(  Oii  <0+1  \  ^  A  max  umin 

l  I  I  J  7  i  = 

\Q!i-l  <0+2  /  OLi 


(23) 


where  e  is  a  small  positive  number  (taken  as  10~6  in  the  numerical  experiments),  and 
Umax  and  wmjn  are  the  maximum  and  minimum  values  of  Uj  over  all  grid  points.  We 
can  see  0  <  0  <  1,  and  0*  =  0( Ax2)  in  smooth  regions  and  0  is  close  to  1  near  a 
discontinuity. 

2  and  fi+1/2  are  modifications  of  fluxes  to  the  anti-diffusive  flux  /^^2. 


m  _  J 

f  /i+1/2 

+  minmod 

j  i+1/2  ) 

1  f(1) 

1 

f  J  i+1/2' 

) 

and 

m  _  J 

\  f i+1/2 

+  minmod 

J  i+1/2  ) 

ii1) 

1 

L  J  i+1/2' 

) 

where  b  = 

Ui  Ui—\ 

A 

^  fi- 1/2  _ 

A 


+  fi- 1/2  fi+l/2i  f i+1/2  fi 


i+1/2 


-1/2  f i+1/2’’  f i+1/2  fi 


-  A 


'i- 1/2 


i+1/2  J 


be  >  0,  |6|  <  |c 

otherwise, 

(24) 

be  >  0,  \b\  <  | c 

otherwise, 

(25) 


It  has  been  numerically  proved  that,  for  a  linear  problem  (i.e.,  ut  +  aux  =  0)  with  a 
piecewise  constant  initial  condition  with  two  values,  the  linear  scheme  (i.e.,  f i+1/2  =  /(u0 
for  a  >  0)  will  not  have  more  than  two  transition  points  between  two  constant  pieces 
for  the  CFL  condition  aX  =  a  At/ Ax  < 

We  refer  to  [33]  for  more  details  about  anti- diffusive  WENO  scheme  in  two  dimensions 
and  will  not  repeat  them  here. 
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2.2  Reaction  operator 


If  there  is  no  smearing  of  discontinuities  in  the  convection  step,  any  ODE  solver  can 
be  used  as  the  reaction  operator.  However,  all  the  standard  shock-capturing  schemes 
will  produce  a  few  transition  points  in  the  shock  when  solving  the  convection  equation. 
These  transition  points  are  usually  responsible  for  causing  incorrect  numerical  results  in 
the  stiff  case.  Thus  we  cannot  directly  apply  a  standard  ODE  solver  at  these  transition 
points. 

Here  we  use  Harten’s  subcell  resolution  technique  in  the  reaction  step.  The  general 
idea  is  as  follows.  If  a  point  is  considered  a  transition  point  of  the  shock,  information 
from  its  neighboring  points  which  are  deemed  not  transition  points  will  be  used  instead. 

The  procedure  can  be  summarized  in  the  following  steps: 

(1)  Use  a  “shock  indicator”  to  identify  cells  in  which  discontinuities  are  believed  to 
be  situated.  We  consider  the  following  minmod-based  shock  indicator  in  [16,  29].  Let 


Si  =  minmod{rq+i  -  rq,  -  rq_ i}, 


(26) 


define  the  cell  It  as  troubled  if  |sj|  >  |si-i|  and  |s*|  >  |sj+i|,  with  at  least  one  being 
a  strict  inequality.  Notice  that  this  troubled  cell-identifying  method  will  only  End  the 
“worst”  cell  inside  a  shock  transition.  That  is,  if  there  are  several  consecutive  transition 
cells,  only  the  worst  one  will  be  identified  as  a  troubled  cell. 

(2)  In  a  troubled  cell  identified  above,  we  continue  to  identify  its  neighboring  cells. 
For  example,  we  can  define  Ii+ 1  as  troubled  if  st+i  \  >  |sj_i|  and  |sj+i|  >  | 1  and 
similarly  define  It-\  as  troubled  if  |s,_i|  >  |sj_2|  and  |sj_i|  >  |si+i|.  If  the  cell  Ii_s  and 
the  cell  Ii+r  ( s,r  >  0)  are  the  first  good  cells  from  the  left  and  the  right  (i.e. ,  J;_s+ 1  and 
Ii+r- 1  are  still  troubled  cells),  we  compute  the  ENO  interpolation  polynomial  p.i_s(x) 
and  Pi+r(x )  for  the  cells  U_.s  and  I^+r^  respectively.  Because  of  the  anti- diffusive  corrector 
in  the  convection  step,  r  and  s  will  not  be  larger  than  2  in  general.  The  modified  cell 
point  value  tq  is  computed  by 


u,  = 


Pi-S(xi),  9  >Xi 

Pi+r{Xi),  9  <  Xi 


(27) 


where  the  location  9  is  determined  by  conservation 


Pi-s(x)dx  + 


ri+ 1/2 


Pi+r(x)dx  =  Ui  Ax. 


(28) 


— 1/2 


Linder  certain  conditions,  it  can  be  shown  that  there  is  a  unique  9  satisfying  Eq. 
(28),  which  can  be  solved  using,  for  example,  a  Newton’s  method.  If  there  is  no  solution 
for  9  or  there  are  more  than  one  solution,  we  choose  fq  =  tq_s. 

(3)  Lise  Uj  instead  of  tq  in  the  ODE  solver  if  the  cell  I,  is  a  troubled  cell. 


For  simplicity,  consider  the  Euler  forward  method 


Eq.  (29)  is  modified  to 


< +1  =  K  +  A  *S(<). 

(29) 

<+I  =  <  +  AtS(Si), 

(30) 

if  the  cell  It  is  a  troubled  cell. 

Here  we  would  like  to  remark  that,  implicit  methods  cannot  be  used  in  this  step 
because  the  troubled  values  u™  need  to  be  modified  explicitly.  However,  there  is  no  small 
time  step  restriction  in  the  explicit  method  used  here,  because  once  the  stiff  points  have 
been  modified,  the  modified  source  term  S(v,j )  is  no  longer  stiff.  Therefore,  a  regular 
CFL  number  is  allowed  in  the  explicit  method. 

For  the  scalar  case,  the  second-order  linearized  trapezoidal  method  is  used  in  the 
numerical  examples 


u?+1 


=  uf  + 


A  tS(ur- 


1  -  f  S' (up 


(31) 


Remark  2.1.  In  general,  the  anti- diffusive  WENO  scheme  can  capture  the  discontinuity 
sharply  with,  at  most,  two  transition  points  inside  the  discontinuity.  Thus  it  does  help 
for  the  stiff  source  term  problems.  For  example,  in  the  numerical  computation,  we  use 
two  immediate  neighbor  cells  (s  —  r  —  1)  for  the  subcell  resolution  procedure.  This 
works  because  the  anti- diffusive  WENO  scheme  provides  a  very  sharp  shock  front  in  the 
convection  step.  However,  if  the  standard  WEN05  is  used  instead,  more  neighbor  cells 
need  to  be  identified  (s,r  >  1)  for  the  subcell  resolution  procedure. 


Remark  2.2.  If  a  multi-step  ODE  solver  is  applied  in  the  reaction  step,  a  modification 
of  the  transition  points  in  each  step  is  required. 


Remark  2.3.  The  proposed  method  is  valid  for  a  general  f(u)  and  a  general  S(u). 


2.3  Numerical  examples  of  ID  scalar  problems 

In  this  section,  we  test  the  proposed  method  on  three  scalar  problems.  The  proposed 
method  uses  a  fifth-order  WENO-Roe  scheme  (WEN05)  with  the  third-order  TVD 
Runge-Kutta  method  (RK3)  as  the  convection  operator,  and  a  linearized  trapezoidal 
method  (31)  based  on  the  subcell  resolution  (SR)  as  the  reaction  operator.  From  now 
on,  we  use  the  notation  WEN05/SR  for  the  proposed  WENO  scheme.  If  the  fifth- 
order  anti- diffusive  WENO  is  used  in  the  convection  step,  the  notation  anti- diffusive 
WEN05/SR  will  be  used. 
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Table  1:  L i  errors  and  orders  of  accuracy  by  the  anti- diffusive  WEN05/SR  at  t  =  0.3 
with  CFL=0.6. 


N 

error 

order 

10 

1.02E-02 

- 

20 

4.06E-04 

4.65 

40 

1.21E-05 

5.07 

80 

3.71E-07 

5.03 

160 

1.18E-08 

4.98 

Example  2.1.  Accuracy  test. 

We  first  test  the  convergence  order  of  the  proposed  anti- diffusive  WEN05/SR  scheme. 
We  consider  Eq.  (4)  with  f{u)  =  u  and  the  source  term 


S(u )  =  —  u  +  sin(27r(a:  —  t )) 


(32) 


and  periodic  boundary  conditions  on  the  computation  domain  x  G  [0,1].  The  exact 
solution  is  u(x,  t )  =  sin(27r(a;  —  t)).  The  errors  and  orders  of  accuracy  are  listed  in  Table 
1.  Since  both  the  Strang  splitting  method  and  also  the  trapezoidal  rule  in  the  reaction 
step  are  second-order  in  time,  we  set  At  =  CFL  x  ( Ax)5// 2  to  achieve  the  fifth  order  in 
space  as  shown  in  Table  1. 

The  next  examples  are  to  show  the  ability  of  the  proposed  schemes  to  deal  with  the 
propagating  shocks. 


Example  2.2.  ID  scalar  test  case  of  a  linear  f{u). 

This  example  is  the  model  problem  of  [23].  Consider  Eq.  (4)  with  f(u)  =  u ,  the 
source  term  given  by  Eq.  (5),  and  the  initial  condition: 


u(x,  0) 


1,  x  <  0.3 
0,  x  >  0.3 


For  this  initial  value  problem,  the  exact  solution  is 


u(x,  t ) 


1,  x  <  t  +  0.3 
0,  x  >  t  +  0.3 


(33) 


(34) 


Analytically,  the  source  term  should  be  always  zero.  However,  if  /i  in  the  source  term 
Eq.  (5)  is  very  large,  the  numerical  errors  of  u  in  the  transition  region  can  result  in  large 
erroneous  values  of  S(u),  which  must  be  corrected. 

We  compare  the  numerical  results  by  the  anti- diffusive  WEN05/SR  and  the  WEN05 
with  splitting  (denoted  by  splitting  WEN05)  in  Figs.  1  and  2,  respectively.  For  each 
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Figure  1:  Results  of  Example  2.2  by  the  anti- diffusive  WEN05/SR  with  N  =  50  at 
t  =  0.3.  Solid  line:  exact  solution;  dashed  line  with  symbols:  computed  solution.  Left: 
H  =  10,  CFL=0.5;  middle:  /i  =  800,  CFL=0.5;  right:  /i  =  10000,  CFL=0.2. 

scheme,  we  test  for  the  cases  of  /x  =  10,  /x  =  800  and  /i  =  10000  with  the  same  mesh 
N  =  50  at  a  final  time  t  =  0.3.  In  the  case  of  fi  =  10,  both  schemes  can  capture 
the  discontinuity  at  the  correct  position  (see  the  left  subplots  of  Figs.  1  and  2).  For 
the  stiffer  case  where  /./  =  800,  the  propagation  speed  of  the  discontinuity  computed 
by  splitting  WEN05  is  qualitatively  slower  than  the  analytical  value  as  shown  by  the 
middle  subplot  of  Fig.  2,  whereas  at  /i  =  10000,  the  discontinuity  solved  by  splitting 
WEN05  does  not  move  at  all.  If  the  mesh  is  sufficiently  refined,  the  splitting  WEN05 
can  capture  the  correct  solution.  However,  for  this  example  in  the  case  where  /i  =  10000, 
at  least  N  =  3000  points  are  needed. 

We  also  note  that  the  anti- diffusive  WENO/SR  is  able  to  produce  correct  results 
with  a  standard  CFL  number.  Even  in  the  very  stiff  case  /i  =  10000,  CFL=0.2  can  be 
used.  But  the  splitting  WEN05  needs  a  very  small  CFL  number  (e.g.,  CFL=0.05)  for 
stability. 

Remark  2.4.  In  this  example,  we  also  show  the  results  by  the  standard  WEN05  without 
splitting  in  Fig.  3.  It  produces  similar  spurious  waves  as  the  splitting  WEN05  with  the 
same  mesh  size,  but  it  requires  a  smaller  CFL  number.  From  now  on  the  numerical 
results  by  the  proposed  scheme  are  only  compared  with  the  results  by  splitting  WEN05. 
The  results  by  standard  WEN05  without  splitting  will  not  be  shown  in  all  of  the  following 
examples. 


Example  2.3.  ID  scalar  test  case  of  a  nonlinear  f(u). 

Consider  a  nonlinear  problem  (4)  with  f{u)  =  \  +  u  and  S(u)  in  (5).  The  initial 
condition  is  (6).  The  discontinuity  has  a  speed  of  |. 
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Figure  2:  Results  of  Example  2.2  by  the  splitting  WEN05  scheme  with  N  =  50  at 
t  =  0.3.  Solid  line:  exact  solution;  dashed  line  with  symbols:  computed  solution.  Left: 
fi  =  10,  CFL=0.5;  middle:  /./  =  800,  CFL=0.1;  right:  /i  =  10000,  CFL=0.05. 


Figure  3:  Results  of  Example  2.2  by  the  standard  WEN05  without  splitting  scheme 
with  IV  =  50  at  t  =  0.3.  Solid  line:  exact  solution;  dashed  line  with  symbols:  computed 
solution.  Left:  /.i  =  10,  CFL=0.5;  middle:  /i  =  800,  CFL=0.1;  right:  /i  =  10000, 
CFL=0.02. 
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Figure  4:  Results  of  Example  2.3  by  different  WENO  schemes  for  /j  =  1000  with  N  =  50 
at  t  =  0.2.  Solid  line:  exact  solution;  dashed  line  with  symbols:  left:  anti-diffusive 
WEN05/SR  with  CFL=0.2;  right:  splitting  WEN05  with  CFL=0.2. 

We  run  the  numerical  experiment  to  t  =  0.2  so  that  the  exact  solution  is  the  same 
as  in  Example  1.  The  discontinuity  moves  to  x  —  0.6  at  t  =  0.2.  We  plot  the  results 
obtained  by  the  anti- diffusive  WEN05/SR  and  splitting  WEN05  schemes  in  the  left 
subplot  and  right  subplot  of  Fig.  4.  Again  N  =  50  and  CFL=0.2  are  used.  Only  the 
anti- diffusive  WEN05/SR  gives  the  correct  numerical  result  (see  the  left  subplot  of  Fig. 
4)  whereas  the  splitting  WEN05  scheme  fails  to  produce  the  correct  shock  speed  in  the 
underresolved  mesh. 

The  result  with  the  stiff  cubic  nonlinear  source  term  (5)  studied  by  LeVeque  and  Yee 
can  be  found  in  [22,  21].  The  cubic  nonlinearity  case  is  slightly  more  complicated  as  one 
or  more  standing  wave  numerical  solutions  can  be  obtained,  depending  on  the  sign  of  fi 
and  the  shock-capturing  method  [22,  21]. 


3  Extension  to  2D  scalar  problems 


It  is  straightforward  to  extend  the  proposed  method  to  the  two-dimensional  scalar  case. 
Consider  the  two-dimensional  scalar  hyperbolic  conservation  law  with  stiff  reaction  term 


u±  +  f(u)x  +  g(u)y  =  S(u), 


(35) 


where  S(u)  is  the  same  as  (5),  i.e.,  S(u)  =  —fm(u  —  a)  (u  —  1),  (0  <  a  <  1),  with 
piecewise  constant  initial  condition 


u(x,y,  0) 


1,  (x,y)  E  fl0  C  M2, 

0,  (x,y)  e  M2  \  fi0, 


(36) 


where  flo  is  a  given  domain  in  M2. 
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We  again  use  the  splitting  method 

un+ 1  =  A  (jpj  R(At)A 

with  the  anti- diffusive  WEN05  as  the  convection  operator  and  the  linearized  trapezoidal 
method  (31)  as  the  ODE  solver  in  the  reaction  step. 

Let  Ii:j  =  [xj_i,xi+i]  x  [yj_i,yj+i\,i  =  1  ,...,Nx,j  =  1,..  .,Ny  be  a  uniform 
partition  of  the  two-dimensional  computational  domain,  with  the  grid  points  x,  = 
\(xi_ i  +  xi+ 1)  and  yj  =  +  Uj+ 1).  Let  u”  be  the  approximated  solution  at 

(xj,  i/j,  tn),  i  =  1, . . . ,  Nx,  j  =  1, . . . ,  iVy.  We  apply  the  subcell  resolution  procedure  in 
the  reaction  step  dimension  by  dimension  in  the  two-dimensional  case. 

(1)  Identify  the  transition  points  by  the  shock  indicator  in  both  x-  and  //-directions. 

Define  the  cell  Il3  as  troubled  in  the  x-direction  if  |s?-|  >  |sf_x  -|  and  |s?-|  >  |sf+1  -| 
with  at  least  one  strict  inequality  where 


Sij  =  minmod{wi+ij  -  ui:j,  ul3  -  ij}.  (38) 

Similarly  we  can  define  the  cell  IVJ  as  troubled  in  the  //-direction  if  141  -  lsij-il and 
I4|  >  1 4.7+1 1  with  at  least  one  strict  inequality  where 

4  =  minmod{wij+i  -  Uij,  ui3  -  ui:j- 1}.  (39) 

If  Iij  is  only  troubled  in  one  direction,  we  apply  the  subccll  resolution  along  this  direc¬ 
tion.  If  1^  is  troubled  in  both  directions,  we  choose  the  direction  which  has  larger  jumps. 
Namely,  if  sf3 1  >  s’l  | ,  subcell  resolution  is  applied  along  the  x-direction,  otherwise  it  is 
done  along  the  //-direction. 

In  the  steps  (2)-(3),  we  assume  take  x-direction,  //-direction  is  similar. 

(2)  Modify  the  point  value  Uij  in  the  troubled  cell  Ii3  by 


f  Pi-sj(xi),  9  >Xi 

(  Pt  -  r.  j  A‘> )  ■  9  Xi 


(40) 


where  the  location  9  is  determined  by  conservation 


Pi-S,j(x)dx  + 


pi+rJ(x)dx 


Ax. 


(41) 


(3)  Use  instead  of  Uij  in  the  ODE  solver  if  the  cell  It]  is  a  troubled  cell,  i.e., 


un+1  =  un ■  _ AtS(u 

av  V  w  At 

-L  o 


ij> 


S'[u 


ijj 


(42) 
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Figure  5:  Results  of  Example  3.1  by  the  anti-diffusive  WEN05/SR  scheme:  y  =  104, 
t  =  0.2,  N  =  50,  CFL=0.1.  Left:  2D  contour;  right:  ID  cross  section  at  y  —  0.5. 


3.1  Numerical  examples  of  2D  scalar  problems 


Example  3.1.  2D  scalar  test  case  of  linear  fluxes. 

Consider  Eq.  (35)  with  f(u)  =  g(u)  =  u  on  the  domain  [0,  l]2,  the  initial  condition 
is  taken  as 


u(x,y,  0) 


1,  x  +  y<  1, 
0,  x  +  y  >  1. 


(43) 


Initially  the  discontinuity  is  located  on  the  diagonal  of  the  square  domain.  It  moves 
at  a  speed  \/2  in  45°  direction.  The  segments  x  =  0  and  y  =  0  are  inflow  boundaries 
condition  and  the  other  two  sides  are  outflow  boundary  condition. 

A  stiff  case  y  =  104  at  time  t  —  0.2  is  considered.  The  left  subplot  of  Fig.  5  shows 
2D  contour  of  the  anti-diffusive  WEN05/SR  with  a  coarse  50  x  50  mesh  and  CFL=0.1. 
The  right  subplot  shows  a  comparison  between  the  anti- diffusive  WEN05/SR  with  the 
splitting  WEN05  using  the  same  mesh  at  y  —  0.5.  At  t  —  0.2,  the  discontinuity  moves  to 
x  =  0.9.  The  proposed  method  is  able  to  capture  the  correct  location  of  the  discontinuity 
with  a  very  coarse  mesh.  However,  the  splitting  WEN05  scheme  fails  to  produce  the 
correct  shock  speed  on  this  underresolved  mesh. 


Example  3.2.  2D  scalar  test  case  of  nonlinear  fluxes. 

In  the  second  example,  we  consider  a  nonlinear  problem  Eq.  (35)  with  f(u)  =  g(u)  = 
Y  on  the  domain  [0,  l]2  with  the  same  initial  condition  (43).  The  boundary  conditions 
are  the  same  as  in  Example  3.1.  In  this  example,  the  discontinuity  moves  at  a  speed 
l/y/2.  A  more  stiff  case  y  =  105  is  tested.  The  left  subplot  of  Fig.  6  shows  the  2D 
contour  results  by  the  anti-diffusive  WEN05/SR  scheme  with  a  coarse  50  x  50  mesh  and 
CFL=0.1.  The  right  subplot  shows  the  ID  cross  section  at  y  —  0.5  by  the  anti- diffusive 
WEN05/SR  and  the  splitting  WEN05  with  the  exact  solution  at  t  =  0.2  where  the 
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Figure  6:  Results  of  Example  2.3  by  the  anti-diffusive  WEN05/SR  scheme:  p  =  105, 
t  =  0.2,  N  =  50,  CFL=0.1.  Left:  2D  contour;  right:  ID  cross  section  at  y  —  0.5. 

discontinuity  moves  to  x  —  0.7.  The  numerical  solution  by  the  proposed  scheme  has 
very  good  agreement  with  the  exact  solution,  but  the  discontinuity  solved  with  the 
splitting  WEN05  does  not  move. 

4  Extension  to  ID  reactive  Euler  equations 

In  this  section,  we  extend  our  approach  to  the  reactive  Euler  equations.  Consider  the 
simplest  one-dimensional  reactive  Euler  equation  with  only  two  chemical  states:  burnt 
gas  and  unburnt  gas.  The  unburnt  gas  is  converted  to  burnt  gas  via  a  single  irreversible 
reaction.  Without  heat  conduction  and  viscosity,  the  system  can  be  written  as 


Pt  +  ( pu)x  =  0, 


(44) 

(45) 

(46) 

(47) 


(pu)t  +  (pu2  +  p)x  =  0, 
Et  +  (u(E  +  p))x  —  0, 


(pz)t  +  (puz)x  =  - K(T)pz , 


where  p  is  the  mixture  density,  u  is  the  mixture  velocity,  E  is  the  mixture  total  energy 
per  unit  volume,  p  is  the  pressure,  z  is  the  mass  fraction  of  the  unburnt  gas,  K  is  the 
chemical  reaction  rate  and  T  is  the  temperature.  The  pressure  is  given  by 


(48) 


where  go  is  the  chemical  heat  released  in  the  reaction.  The  temperature  is  defined  as 


(49) 


P 
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The  reaction  rate  K(T )  is  modeled  by  an  Arrhenius  law 


K(T)  =  Kaex (50) 

where  A'0  is  the  reaction  rate  constant  and  Tign  is  the  ignition  temperature.  The  reaction 
rate  may  be  also  modeled  in  the  Heaviside  form 

K<r>-{l,tTT*<£’  <51> 

where  e  is  the  reaction  time  and  1/e  is  roughly  equal  to  Kq. 

We  treat  (44)-(47)  similarly  as  the  scalar  case. 

4.1  Convection  operator 

In  the  scalar  problem  Eq.  (4),  we  have  applied  the  anti-diffusive  WEN05  scheme  as  the 
convection  operator  because  the  discontinuous  wave  is  a  contact  discontinuity.  However, 
the  Chapman- Jouguet  (C-J)  detonation  wave  is  a  shock  followed  by  a  reaction.  In  gen¬ 
eral,  it  is  not  safe  to  apply  the  anti- diffusive  technique  to  a  shock,  since  it  may  generate 
an  entropy-violating  solution.  Therefore,  we  do  not  apply  the  anti- diffusive  sharpening 
procedure  here.  This  is  not  a  problem  because  the  standard  WEN05  scheme  is  already 
able  to  capture  the  shock  very  sharply  (better  than  its  capability  to  capture  contact 
discontinuities).  In  the  system  case,  we  use  WEN05  with  Lax-Friedrichs  flux  splitting 
(WENO-LF)  and  the  local  characteristic  decomposition  with  RK3  in  time  discretization 
as  the  convection  operator  in  the  reactive  Euler  problems.  We  refer  to  [19]  for  more 
details. 

4.2  Reaction  operator 

The  reaction  step  for  the  system  case  is  slightly  different  from  the  scalar  case  because 
there  are  more  component  variables  (pz  and  T )  involved  in  the  source  term.  The  key 
point  here  is  to  identify  transition  points  correctly  and  to  extrapolate  the  variables  pz 
and  T. 

(1)  To  apply  step  (1)  in  Section  2.2,  we  need  to  choose  one  variable  to  examine.  Note 
that  in  a  detonation  wave,  the  pressure,  temperature,  and  density  all  have  a  reaction 
zone  (like  an  “overshoot”)  and  a  shock  zone.  Only  the  mass  fraction  z  has  a  clean  single 
shock  wave.  This  can  be  seen  from  the  mass  fraction  equation.  Eliminating  the  density 
from  Eq.  (47)  by  using  Eq.  (44),  we  obtain 

zt  +  uzx  =  - K(T)z ,  (52) 
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which  is  of  the  scalar  type  Eq.  (4).  This  helps  us  identify  transition  points  by  the 
variable  z.  Define  the  cell  E  as  troubled  if  |sj|  >  |sj_i|  and  |s*|  >  |sj+i|  (with  at  least 
one  strict  inequality)  where 

Si  =  minmod{zj+i  -  zh  Zi  -  J.  (53) 


(2)  After  a  troubled  cell  I,  is  identified,  first  find  the  shock  location  6  by  solving  the 
conservation  Eq.  (28)  with  the  variable  u  taken  as  the  total  energy  E 

rxi+ 1/2 

Pi-S(x]  E)dx  +  /  pi+r(x;  E)dx  =  Ei  Ax,  (54) 

1/2  JO 

where  the  ENO  interpolation  polynomials  pi(x’,E )  are  computed  based  on  values  of  E. 
The  energy  E  is  chosen  because  it  is  a  conserved  variable.  We  assume  the  shock  locations 
are  the  same  for  all  variables.  Then  we  extrapolate  the  variables  p,  z  and  T  separately. 
The  new  mass  fraction  z,  temperature  T  and  density  p  are  obtained  from  the  ENO 
interpolation  polynomials 

f  ^  =  Pi-a{xi ;  z),  Tj=  Pisixp  T ),  pi=  Pi-S(xi\ p),  if  6  >  xt 
\  Zi=  Pi+r(Xi ;  z),  Ti  —  Pi+rixp  T ),  Pi  =  pi+r(xi ;  p),  if  9  <  Xi 

Remark  4.1.  s  —  r  —  1  works  well  in  all  the  numerical  examples  for  the  system  case. 


Remark  4.2.  Observe  that  the  mass  fraction  z  has  values  0  or  1  for  the  burnt  gas  arid 
unburnt  gas,  respectively.  Values  between  0  and  1  denote  the  gas  changing  from  unburnt 
to  burnt.  However,  in  stiff  reaction  problems,  the  reaction  is  very  fast  and  the  reaction 
zone  is  much  smaller  than  the  grid  size  for  an  underresolved  mesh.  Thus  in  stiff  reaction 
problems,  the  grid  values  of  z  can  be  simplified  to  have  only  two  values  0  and  1,  but  no 
middle  values,  i.e., 

0,  9  >  Xi 

1,  9  <  Xi 

instead  of  the  values  obtained  from  the  ENO  polynomial  extrapolation. 


(3)  For  simplicity,  we  use  the  explicit  Euler  method  as  the  ODE  solver  in  the  reaction 
step 

(pz)^'  =  (pz^  +  AtS(fl,pl,zi).  (57) 

In  general,  a  regular  CFL=0.1  can  be  used  in  the  proposed  scheme  to  produce  a 
stable  solution.  But  the  solution  is  very  coarse  in  the  reaction  zone  because  of  the 
underresolved  mesh  in  time.  In  order  to  obtain  more  accurate  results  in  the  reaction 
zone,  we  evolve  one  reaction  step  via  Nr  sub  steps,  i.e., 


in  some  numerical  examples. 
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4.3  Numerical  examples  of  one-dimensional  detonation  waves 


In  this  section,  we  test  the  proposed  method  on  five  examples  of  one- dimensional  detona¬ 
tion  waves.  The  first  example  uses  the  Arrhenius  law  (50).  The  next  four  examples  are 
based  on  the  Heaviside  model.  The  proposed  method  uses  a  fifth-order  WENO-LF  with 
RK3  as  the  convection  operator,  and  an  explicit  Euler  based  on  the  subcell  resolution 
as  the  reaction  operator. 


Example  4.1.  C-J  detonation  wave  (Arrhenius  case). 

The  first  example  uses  the  Arrhenius  source  term  (50)  and  has  been  studied  in 
[17,  32],  The  initial  values  consist  of  totally  burnt  gas  on  the  left-hand  side  and  totally 
unburnt  gas  on  the  right-hand  side.  The  density,  velocity,  and  pressure  of  the  unburnt 
gas  are  given  by  pu  —  1,  uu  —  0  and  pu  =  1.  The  heat  release  q0  =  25  and  the  ratio 
of  specific  heats  is  set  to  7  =  1.4.  We  consider  the  ignition  temperature  Tign  =  25  and 
Kq  =  164180.  We  can  obtain  the  C-J  initial  state  for  the  unburnt  gas  by,  for  example, 
[12] 


Pb  = 

-b+  (b2  -  c)1/2, 

(59) 

pb  = 

Pu[Pb( 7  +  1)  ~Pu] 

(60) 

IPb 

scj  = 

[pu'U'u  (SIPbpb)  ^  ]/ Pm 

(61) 

ub  = 

SCJ  -  ( IPb/pb)11 2, 

(62) 

where 


b  =  ~Pu  ~  Puqoh  -  1), 

c  =  Pu  +  2(7  -  l)pup«9o/(7  +  !), 


and  scj  is  the  speed  of  the  C-J  detonation  wave.  In  this  example,  Scj  =  7.1247. 

The  computational  domain  is  [0,  30].  Initially,  the  discontinuity  is  located  at  x  =  10. 
At  time  t  =  1.8,  the  detonation  wave  has  moved  to  x  —  22.8.  The  reference  solution  is 
computed  by  the  standard  WEN05  scheme  with  N  =  10000  (Ax  =  0.003),  CFL=0.05. 

Figures  7-10  show  the  pressure,  temperature,  density  and  mass  fraction  comparison 
results  between  WEN05/SR  method  with  the  splitting  WEN05  method.  Only  N  =  50 
(Ax  =  0.6)  and  CFL=0.1  are  used  in  WEN05/SR.  Clearly,  our  WEN05/SR  method  is 
able  to  capture  the  correct  propagation  speed  of  the  detonation  wave  with  this  coarse 
mesh,  while  the  splitting  WEN05  with  a  much  finer  mesh  N  =  300  produces  spurious 
numerical  results.  There  are  small  downstream  wiggles  located  around  x  —  8  in  Figs. 
Figures  7-9  which  are  numerical  artifacts.  They  become  smaller  as  the  mesh  refines  and 
will  move  away  after  a  period  of  time. 

We  remark  that  our  method  can  use  fewer  points  than  the  previous  methods  in  [17] 
and  [32]  to  obtain  similar  results.  The  reason  may  be  due  to  the  high  order  accuracy  of 
the  spatial  scheme  in  the  convection  step. 
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Figure  7:  Computed  pressure  for  the  Arrhenius  Example  4.1  at  t  —  1.8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1;  right:  splitting  WEN 05  with  N  =  300,  CFL=0.1. 
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Figure  8:  Computed  temperature  for  the  Arrhenius  Example  4.1  at  t  —  1.8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1;  right:  splitting  WEN05  with  N  =  300,  CFL=0.1. 
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Figure  9:  Computed  density  for  the  Arrhenius  Example  4.1  at  t  —  1.8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1;  right:  splitting  WEN 05  with  N  =  300,  CFL=0.1. 


Figure  10:  Computed  mass  fraction  for  the  Arrhenius  Example  4.1  at  t  —  1.8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1;  right:  splitting  WEN05  with  N  =  300,  CFL=0.1. 
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Example  4.2.  C-J  detonation  wave  (Heaviside  case). 

In  this  example  the  chemical  reaction  is  modeled  by  the  Heaviside  form.  This  example 
is  taken  from  [13,  5,  2], 

Consider  the  following  parameter  values  in  CGS  units: 

7  =  1.4,  q0  =  0.5196  x  1010,  ^  =  0.5825  x  1010,  Tign  =  0.1155  x  1010. 

The  computational  domain  is  [0,0.05].  The  initial  conditions  are  given  by 


f  (pb,ub,pb,  0)  x  <  0.005 

\  (1.201  x  10”3,  0,8.321  x  105,1)  x  >  0.005  ’ 


(63) 


where  pb,  ub  and  pb  are  computed  by  Eqs.  (59)-(62).  From  Eq.  (61),  the  speed  of  the 
detonation  wave  in  this  example  is  Dqj  =  1.088  x  105.  In  this  example,  the  width  of 
the  reaction  zone  is  approximately  5  x  10”5  (see  [5]  and  [13]). 

The  reference  solution  is  computed  by  the  standard  WEN05  scheme  with  N  =  5000 
points  ( Ax  =  10”5)  and  CFL=0.05.  The  solutions  are  run  to  time  t  —  3  x  10”'.  The 
wave  moves  to  x  —  0.03764.  The  pressure,  temperature,  density  and  mass  fraction 
results  are  plotted  in  Figs.  11-14.  IV  =  50  and  CFL=0.1  are  used  in  WENO/SR.  In 
this  example,  10  sub  reaction  steps  are  involved  in  each  time  step  evolution  in  order  to 
produce  more  accurate  results  in  the  reaction  zone. 

Again  we  compare  the  results  by  WEN05/SR  (left  subplots)  with  splitting  WEN05 
(right  subplots).  We  can  see  WEN05/SR  with  IV  =  50  is  able  to  capture  the  correct 
detonation  speed.  However,  splitting  WEN05  with  N  =  300  still  produces  wrong  nu¬ 
merical  results  no  matter  how  small  the  time  step  is  (the  results  with  smaller  time  steps 
are  not  shown  here  to  save  space). 

Example  4.3.  A  strong  detonation  (Heaviside  case). 

Here  is  another  detonation  problem  which  is  also  from  [2] .  The  computational  domain 
is  [0,0.05].  The  initial  data  are 


(  (pi,ut,pi,  0)  x  <  0.005 
1  (pr,  Ur,pr,  1)  X  >  0.005 


(64) 


where  ui  =  9.162  x  104  >  ub,  p\  =  pb,  Pi  =  Pb,  and  the  right  state  is  the  same  as  in 
Example  4.2.  The  exact  solution  contains  a  right-moving  strong  detonation,  a  right- 
moving  contact  discontinuity  and  a  stationary  shock. 

The  reference  solutions  are  computed  by  standard  WEN05  with  N  =  5000  (Ax  = 
1  x  10”5)  and  CFL=0.05.  We  display  the  numerical  results  by  WEN05/SR  with  N  =  50, 
CFL=0.1  and  Nr  =  10  at  t  =  2  x  10”'.  The  pressure,  temperature,  density  and  mass 
fraction  results  are  plotted  in  the  left  subplots  of  Figs.  15-18.  We  also  compute  the 
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Figure  11:  Computed  pressure  for  the  Heaviside  Example  4.2  at  t  —  3  x  10  ' .  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  300,  CFL=0.01. 


Figure  12:  Computed  temperature  for  the  Heaviside  Example  4.2  at  t  =  3  x  10~'. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  300, 
CFL=0.01. 
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Figure  13:  Computed  density  for  the  Heaviside  Example  4.2  at  t  —  3  x  10”'.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR. 
with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  300,  CFL=0.01. 


Figure  14:  Computed  mass  fraction  for  the  Heaviside  Example  4.2  at  t  =  3  x  10”'. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  300, 
CFL=0.01. 
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Figure  15:  Computed  pressure  for  the  Heaviside  Example  4.3  at  t  —  2  x  10  ' .  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  50,  CFL=0.01. 

results  by  the  splitting  WEN05  with  N  —  50  and  CFL=0.01  in  the  right  subplots  of 
Figs.  15-18.  We  can  see  WEN05/SR  is  able  to  capture  the  correct  shock  speed  and 
other  waves  in  a  very  coarse  mesh.  But  the  splitting  WEN05  with  the  same  mesh 
produce  spurious  waves. 

Example  4.4.  Collision  of  a  detonation  with  a  rarefaction  wave  (Heaviside 
case). 

Next,  we  consider  a  one- dimensional  detonation  problem  involving  a  collision  with  a 
rarefaction  wave.  This  example  is  taken  from  [3]  and  [18].  The  parameters  are  taken  as 
7  =  1.2,  q0  =  50,  i  =  230.75  and  Tign  =  3. 

The  computational  domain  is  [0,100].  The  initial  data  are 

{(pi,uhpi,  0)  x  <  10 

(Pm,Um,Pm,  0)  10  <  X  <  20  (65) 

(pr,  Ur,pr,  1)  X  >  20 

where  pi  =  2,  ui  =  4,  pi  =  40,  pm  =  3.64282,  pm  =  54.8244,  um  =  6.2489,  pr  —  1,  ui  —  0 
and  pi  =  1. 

The  exact  solution  contains  a  right-moving  strong  detonation,  a  right  moving  rar¬ 
efaction  wave,  a  right  moving  contact  discontinuity,  and  a  left  moving  rarefaction  wave. 
After  some  time,  the  right  moving  rarefaction  wave  will  catch  up  with  the  detonation 
wave.  We  consider  the  solutions  before  the  collision  of  the  detonation  with  the  rarefac¬ 
tion  at  t  —  2  and  the  solutions  after  the  collision  at  t  —  8.  The  reference  solutions  are 
computed  by  standard  WEN05  with  N  =  10,000  (Ax  =  0.01)  and  CFL=0.3. 

Figure  19  shows  pressure,  temperature,  density  and  mass  fraction  results  by  WEN05/SR 
with  N  =  100  and  CFL=0.1  at  t  —  2.  Before  the  collision,  both  the  WEN05/SR  and 
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Figure  16:  Computed  temperature  for  the  Heaviside  Example  4.3  at  t  =  2  x  10  '. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  with  N  =  50, 
CFL=0.01. 
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Figure  17:  Computed  density  for  the  Heaviside  Example  4.3  at  t  =  2  x  10~7.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  WEN05  with  N  =  50, 
CFL=0.01. 


26 


Figure  18:  Computed  mass  fraction  for  the  Heaviside  Example  4.3  at  t  =  2  x  10~7. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  50,  CFL=0.1,  Nr  =  10;  right:  splitting  WEN05  WEN05  with 
N  =  50,  CFL=0.01. 


the  splitting  WEN05  method  can  capture  the  correct  shock  speed  on  the  mesh  with  100 
points  where  the  results  by  the  splitting  WEN05  are  not  shown  here.  However,  after 
the  collision  at  t  —  8  (see  Figs.  20-23),  the  splitting  WEN05  cannot  capture  the  correct 
shock  location  and  produce  spurious  numerical  waves  around  the  detonation  which  can 
be  easily  seen  in  pressure,  temperature  and  mass  fraction  results  (see  right  subplots  of 
Figs.  20,  21,23).  But  the  spurious  waves  are  very  small  in  the  density  result  which 
appear  in  the  bottom  corner  around  the  detonation  (see  right  subplots  of  Figs.  22.  The 
proposed  WEN05/SR  scheme  is  able  to  capture  the  correct  shock  speed  and  other  waves 
in  a  very  coarse  mesh. 

Example  4.5.  A  detonation  interacting  with  an  oscillatory  profile  (Heaviside 
case). 

The  last  one-dimensional  detonation  problem  involves  a  collision  with  an  oscillatory 
profile.  This  example  is  also  taken  from  [3].  The  parameters  7,  g0  and  Tign  are  the  same 
as  Example  4.4  except  -  =  1000. 

The  computational  domain  is  [0 ,  27t]  .  The  initial  data  are 


(p,  u,  p,  z ) 


(pi,ui,pi,  0)  x  <  f 

(pr,Ur,pr,  1)  X  >  f 


(66) 


where  pi  =  1.79463,  ui  =  3.0151,  pi  =  21.53134,  pr  =  1.0  +  0.5  sin  2x,  ui  =  0  and  pi  =  1. 

The  reference  solutions  are  computed  by  splitting  WEN05  with  N  =  10,000  and 
CFL=0.3.  The  numerically  computed  pressure,  temperature,  density  and  mass  fraction 
at  t  —  7t/5  are  plotted  in  Figs.  24-27  separately,  where  the  left  subplots  are  computed 
by  WEN05/SR  with  N  =  200  and  CFL=0.1,  and  the  right  subplots  are  computed  by 
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Figure  19:  Results  for  the  Heaviside  Example  4.4  at  t  —  2.  Solid  line:  reference  solution. 
Dashed  line  with  symbols:  numerical  solution  of  WEN05/SR  with  N  =  100,  CFL=0.1. 
Top  left:  pressure;  top  right:  temperature;  bottom  left:  density;  bottom  right:  mass 
fraction. 


Figure  20:  Computed  pressure  for  the  Heaviside  Example  4.4  at  t  =  8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  100,  CFL=0.1;  right:  splitting  WEN05  with  N  =  100,  CFL=0.1. 
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Figure  21:  Computed  temperature  for  the  Heaviside  Example  4.4  at  t  —  8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  100,  CFL=0.1;  right:  splitting  WEN05  with  N  =  100,  CFL=0.1. 


Figure  22:  Computed  density  for  the  Heaviside  Example  4.4  at  t  =  8.  Solid  line: 
reference  solution.  Dashed  line:  numerical  solution.  Left:  WEN05/SR  with  N  =  100, 
CFL=0.1;  right:  splitting  WEN05  with  N  =  100,  CFL=0.1. 


29 


Figure  23:  Computed  mass  fraction  for  the  Heaviside  Example  4.4  at  t  =  8.  Solid  line: 
reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left:  WEN05/SR 
with  N  =  100,  CFL=0.1;  right:  splitting  WEN05  with  N  =  100,  CFL=0.1. 

splitting  WEN05  with  N  =  200  and  CFL=0.1.  We  can  see  the  proposed  WEN05/SR 
is  able  to  handle  the  interactions  between  the  detonation  and  the  oscillatory  wave  in 
a  very  coarse  mesh,  while  the  splitting  WEN05  scheme  produces  unphysical  solutions 
around  the  detonation  shock. 


5  Extension  to  2D  reactive  Euler  equations 

Next,  we  extend  the  proposed  method  to  the  two-dimensional  reactive  Enler  equations. 
The  considered  two-dimensional  problem  is  the  extension  of  the  one-dimensional  prob¬ 
lem,  again  modeling  the  reaction  with  two  chemical  states  and  one  reaction.  The  gov¬ 
erning  equations  are 


Pt  +  (pu)x  +  ( pv)y 

=  0 

(67) 

(. PU)t  +  ( PU 2  +  P)X  +  ( pUV)y 

=  0 

(68) 

(pv)t  +  ( PUV)X  +  (pV2  +  p)y 

=  0 

(69) 

+  (u(E  +  p))x  +  (v(E  +  p))y 

=  0 

(70) 

G pz)t  +  {puz)x  +  {pvz)y 

~K(T)pz, 

(71) 

where  p(x,y,t )  is  the  mixture  density,  u(x,y,t )  and  v(x,y,t )  are  the  mixture  x-  and  y- 
velocities,  E(x,  y,  t )  is  the  mixture  total  energy  per  unit  volume,  p(x,  y,  t )  is  the  pressure, 
z(x,  y,  t )  is  the  mass  fraction  of  the  unburnt  gas,  K(T )  is  the  chemical  reaction  rate  and 
T(x,y,t )  is  the  temperature.  The  pressure  is  given  by 

P  =  (7  -  1)(£  -  +  v2)  ~  qopz),  (72) 
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Figure  24:  Computed  pressure  results  for  the  Heaviside  Example  4.5  at  t  —  tt/5. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  200,  CFL=0.1;  right:  splitting  WEN05  with  N  —  200,  CFL=0.1. 
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Figure  25:  Computed  temperature  results  for  the  Heaviside  Example  4.5  at  t  =  7t/5. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  iV  =  200,  CFL=0.1;  right:  splitting  WEN05  with  N  —  200,  CFL=0.1. 
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Figure  26:  Computed  density  results  for  the  Heaviside  Example  4.5  at  t  —  7t/5. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  TV  =  200,  CFL=0.1;  right:  splitting  WEN05  with  N  —  200,  CFL=0.1. 
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Figure  27:  Computed  mass  fraction  results  for  the  Heaviside  Example  4.5  at  t  —  7t/5. 
Solid  line:  reference  solution.  Dashed  line  with  symbols:  numerical  solution.  Left: 
WEN05/SR  with  N  =  200,  CFL=0.1;  right:  splitting  WEN05  with  N  —  200,  CFL=0.1. 
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where  the  temperature  T  =  -p  and  qo  is  the  chemical  heat  released  in  the  reaction.  The 
source  term  is  modeled  as  in  the  one-dimensional  case.  For  simplicity,  we  only  consider 
the  Heaviside  source  term  (51). 

In  the  convection  step,  we  use  fifth-order  WENO-LF  with  RK3  time  discretization. 

In  the  reaction  step,  we  apply  the  subcell  resolution  procedure  dimension  by  dimen¬ 
sion. 

(1)  Identify  troubled  cell  1,-j  in  both  x-  and  //-directions  by  applying  the  shock  indi¬ 
cator  to  z. 

Assuming  1^  is  troubled  in  the  x-direction,  we  apply  subcell  resolution  along  the 
x-direction. 

(2)  Modify  the  point  value  %  Tl}  and  ptJ  in  the  troubled  cell  ItJ  by  the  ENO  inter¬ 
polation  polynomials 


f  Zij  =  Pi-a,j{xi ;  z),  fij  =  Pi_sj(xi]  T),  pij  =  Pi_sj(xi]  p ),  if  9  >  x* 

\  Aj  Pi+r.  j  (Xj ,  Z) ,  Tij  Pi-\-r.  j  (X{ ,  T) ,  Pij  Pi+r,j  (Xi ,  p),  if  0  Xj 


where  the  location  6  is  determined  by  the  conservation  of  energy  E 


Pi—S,j  (x;  E)dx  + 


Pi+rtj(x)  E)dx 


EijAx. 


For  simplicity,  in  the  considered  stiff  problem,  the  value  of  can  be  taken  as 


(73) 


(74) 


f  0,  9  >  Xi 

\  1,  9  <  Xi 


(75) 


(3)  For  simplicity,  explicit  Euler  is  used  as  the  ODE  solver. 


5.1  Numerical  examples  of  2D  detonation  waves 

Example  5.1.  A  2D  detonation  wave. 

This  example  is  taken  from  [2],  The  chemical  reaction  is  modeled  by  the  Heaviside 
form  with  the  same  parameters  qo,  -  and  Tign  as  in  Example  4.2.  Consider  a  two- 
dimensional  channel  of  width  0.005,  the  upper  and  lower  boundaries  are  solid  walls. 
The  computational  domain  is  [0,0.025]  x  [0,0.005].  The  initial  conditions  are 


(p,u,v,p,z) 


(pi,ui,0,pi,0),  if  x<$(y), 
(pr,  Ur,  0,pr,  1),  if  X  >  £(p), 


(76) 


where 


£(y)  = 


0.004  \y  -  0.0025|  <  0.001, 

0.005  -  | y-  0.0025|  j y  -  0.0025|  <  0.001, 


(77) 
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and  ui  =  8.162  x  104  >  Ub ,  pi  =  Pb  and  pi  =  pb ■  Ub,  Pb,  Pb  and  the  right  state  are  as  in 
Example  4.2. 

Similar  problems  are  also  computed  in  [14],  One  important  feature  of  this  solution 
is  the  appearance  of  triple  points,  which  travel  in  the  transverse  direction  and  reflect 
from  the  upper  and  lower  walls.  A  discussion  of  the  mechanisms  driving  this  solution  is 
given  in  [20]. 

Figures  28-29  show  density  contours  computed  by  WEN05/SR  with  500  x  100  (Ax  = 
Ay  =  5  x  10  5),  CFL=0.1  and  Nr  =  2  at  eighteen  evolutionary  times  from  t  —  0  to 
t  =  1.7  x  10”7.  We  can  see  the  movement  of  the  triple  points.  The  same  case  by 
WEN05/SR  with  a  much  coarser  grid  200  x  40  (Ax  =  Ay  =  1.25  x  10”4)  with  CFL=0.1 
and  Nr  =  2  at  three  evolutionary  times  are  shown  in  Fig.  30.  We  can  see  WEN05/SR 
with  the  very  coarse  200  x  40  mesh  can  still  capture  the  correct  shock  location,  although 
the  shocks  are  smeared  due  to  the  lack  of  resolution.  It  is  more  apparent  to  compare 
the  computed  results  with  the  reference  solution  in  a  ID  cross  section.  The  reference 
solutions  are  computed  by  standard  WEN05  with  2000  x  400  grid  points  and  CFL=0.3. 
The  results  by  WEN05/SR  and  the  splitting  WEN05  are  compared  with  the  same 
mesh  200  x  40  and  CFL=0.005.  Figures  31-34  show  the  ID  cross  section  at  y  =  0.0025 
at  evolutionary  times  t  =  2  x  10“8,  t  =  6  x  1CT8,  t  —  1.4  x  10“7  and  t  =  1.7  x  1CT7 
separately,  where  the  left  subplots  are  computed  by  WEN05/SR  and  the  right  subplots 
are  by  splitting  WEN05.  We  can  see  WEN05/SR  has  excellent  agreement  with  the 
reference  solutions  except  it  cannot  capture  the  waves  sharply  due  to  the  underresolved 
mesh.  However  the  splitting  WEN05  method  produces  spurious  waves  in  front  of  the 
detonation  shock  at  the  beginning  time  t  =  2  x  10~8  (right  subplot  of  Fig.  31)  and  after 
that  the  solutions  move  at  a  wrong  speed  (right  subplots  of  Figs.  32-34). 


6  Concluding  remarks 

A  new  high  order  finite  difference  scheme  with  subcell  resolution  for  hyperbolic  conser¬ 
vation  laws  with  stiff  source  terms  has  been  developed.  This  method  utilizes  a  fractional 
step  approach  with  the  freedom  in  choosing  any  spatial  high-resolution  shock-capturing 
schemes  and  temporal  discretizations.  In  the  convection  step,  any  spatial  high-resolution 
scheme  can  be  used.  In  the  reaction  step,  any  explicit  ODE  solver  can  be  used  with 
the  transition  points  reconstructed  by  Harten’s  ENO  subcell  resolution.  The  proposed 
method  has  high  order  accuracy  in  space  for  smooth  flows.  It  is  able  to  capture  the  cor¬ 
rect  discontinuity  speed  on  very  coarse  meshes  and  with  a  reasonable  CFL  number.  It 
can  be  used  for  both  stiff  and  non-stiff  problems.  Extensive  numerical  examples  for  one- 
and  two-dimensional  scalar  problems  and  one-  and  two-dimensional  detonation  waves 
demonstrate  the  robustness  of  the  method. 

From  the  numerical  experiments,  further  containment  of  numerical  dissipation  in 


34 


0.007  0.008 
x 


0.008  0.009  0.01 
x 


0.009  0.01  0.011 
x 


Figure  28:  Computed  density  for  Example  5.1:  WEN05/SR  with  500  x  100,  CFL=0.1 
and  Nr  —  2  at  nine  different  evolutionary  times  t  —  0,  t  —  1  x  10~8,  t  —  2  x  10~8, 
t  =  3  x  10"8,  t  =  4x  10~8,  t  =  5  x  10"8,  t  =  6  x  10~8,  t  =  7  x  10"8  and  t  =  8  x  10“8. 
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Figure  29:  Computed  density  for  Example  5.1:  WEN05/SR  with  500  x  100,  CFL=0.1 
and  Nr  —  2  at  nine  different  evolutionary  times  t  =  9  x  10-8,  t  —  lx  10-',  t  —  1.1  x  10-7, 
t  =  1.2  x  10~7,  t  =  1.3  x  10“7,  t  =  1.4  x  10"7,  f  =  1.5  x  10”7,  i  =  1.6  x  10"7  and 
t  =  1.7  x  10~7. 
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Figure  30:  Density  results  of  Example  5.1:  WEN05/SR  with  200  x  40,  CFL=0.1  and 
Nr  =  2  at  t  =  1.5  x  IQ’7,  t  =  1.6  x  10~7  and  t  =  1.7  x  KT7. 


Figure  31:  ID  cross-section  of  Example  5.1  at  t  =  2  x  10~8  by  different  WENO  schemes 
with  200  x  40.  Solid  line:  reference  solution;  dashed  line:  numerical  solution.  Left: 
WEN05/SR  with  CFL=0.1,  Nr  =  2;  right:  splitting  WEN05  with  CFL=0.05. 
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Figure  32:  ID  cross-section  of  Example  5.1  at  t  —  6  x  ICE8  by  different  WENO  schemes 
with  200  x  40.  Solid  line:  reference  solution;  dashed  line:  numerical  solution.  Left: 
WEN05/SR  with  CFL=0.1,  Nr  =  2;  right:  splitting  WEN05  with  CFL=0.05. 


Figure  33:  ID  cross-section  of  Example  5.1  at  t  —  1.4  x  10''  by  different  WENO 
schemes  with  200  x  40.  Solid  line:  reference  solution;  dashed  line:  numerical  solution. 
Left:  WEN05/SR  with  CFL=0.1,  Nr  =  2;  right:  splitting  WEN05  with  CFL=0.05. 
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Figure  34:  ID  cross-section  of  Example  5.1  at  t  —  1.7  x  ICE'  by  different  WENO 
schemes  with  200  x  40.  Solid  line:  reference  solution;  dashed  line:  numerical  solution. 
Left:  WEN05/SR  with  CFL=0.1,  Nr  =  2;  right:  splitting  WEN05  with  CFL=0.05. 


existing  high  order  shock-capturing  schemes  can  defer  the  onset  of  wrong  propagation 
speeds  of  discontinuities  to  certain  degree.  However,  the  need  to  contain  the  spreading 
of  the  discontinuity  front  is  the  key  to  overcoming  the  difficulty.  In  future  work,  we 
will  extend  this  approach  to  more  general  chemical  reaction  models  including  multiple 
reaction  models. 

Due  to  the  fact  the  current  approach  is  only  second  order  in  time  due  to  the  splitting 
method,  we  will  also  consider  developing  a  non-splitting  method  with  an  explicit  RK 
scheme  with  the  subccll  resolution  applied  to  the  source  term.  However  this  is  not  a 
trivial  task.  Straight  application  of  the  subccll  resolution  to  the  source  term  in  each 
inner  stages  of  RK  scheme  does  not  work.  Because  in  the  system  case,  the  source  term 
should  choose  the  correct  value  for  each  time  step.  If  there  are  three  inner  stages  in  RK 
scheme  and  each  stage  takes  a  different  value  for  the  source,  the  convex  combination  of 
them  may  lead  to  wrong  results.  We  will  investigate  this  more  in  the  future. 
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